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The transition of the flow in a duct of square cross-section is studied. Like in the similar 
case of the pipe flow, the motion is linearly stable for all Reynolds numbers; this flow is 
thus a good candidate to investigate the 'bypass' path to turbulence. Initially the so-called 
'linear optimal perturbation problem' is formulated and solved, yielding optimal distur- 
bances in the form of longitudinal vortices. Such optimals, however, fail to elicit a signifi- 
cant response from the system in the nonlinear regime. Thus, streamwise-inhomogeneous, 
sub-optimal disturbances are focussed upon; nonlinear quadratic interactions are imme- 
diately evoked by such initial perturbations and an unstable streamwise-homogeneous 
large amplitude mode rapidly emerges. The subsequent evolution of the flow, at a value 
of the Reynolds number at the edge between fully developed turbulence and relaminariza- 
tion, shows the alternance of patterns with two pairs of large scale vortices near opposing 
parallel walls. Such edge states bear a resemblance to optimal disturbances. 
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1. Introduction 

Transition to turbulence in ducts is still an unsolved issue despite the 120-plus years 
since the observations by Osborne Reynolds that led to the definition of a similarity 
parameter, the ratio of the viscous to the convective time scale, capable of broadly sepa- 
rating the cases where the flow state was laminar from those where turbulence prevailed. 
Recent years have seen a resurgence of interest in the topic, spurred by new developments 
in linear and nonlinear stability theories. As is now well known, classical small perturba- 
tion t heory is unc a pable to provide an e x planation for the onset of tran sition in ducts and 
pipes dcii (|l965l) : ISalwen ei all (|l980l ): iTatsumi fc Yoshimural (|l990f )). Current under- 
standing ascribes the failure of classical theory to its focus on the asymptotic behavior 
of individual modes; when a small disturbance composed by a weighted combination 
of linear eigenfunctions is considered, there is the potential for very large short-time 
amplification of perturbation energy, even in no minally stable flow conditions. This be- 
havior has been reported by iLandahll (|1980 ) and iBoberg fc Brosal (|l988) and has been 
given the name of algebraic instability (and later 'transient growth theory'), since the 
initial rapid growth in time of small disturbances goes like t. The property is related to 
the non-normality of the linearized stability operator (which does not commute with its 
adjoint). 

Despite the appeal and elegance of transient growth theory, it was realized that the fully 
nonlinear Navier-Stokes equations need to be u s ed to under s tand transition phenomena . 
In this context, we mention th e work oflNagatal (ll990lll997l) .IWaleffel (Il997lll998il2003h. 
Faisst fc Eckhardtl (|2003l ) and IWedin fc Kerswelll (|2Q04| ). These works present traveling 
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wave and equilibrium solutions of the Navier-Stokes equation for channel and pipe flows 
that a re possibly related to t ransition. Experimental investigations along these lines are 
due to lHof e<aLl(|2004Ll2005l) . 



In dynamical systems' terminology it is argued that unstable travelling waves appear 
through saddle node bifurcations in phase space; the travelling waves found so far are all 
saddle points with low dimensional unstable manifolds. The saddles act by attracting the 
flow from the vicinity of the laminar state, and then repelling it away. For transitional 
or turbulent, yet moderate, values of Re the flow wanders in phase space between few 
repelling states, spending much time in their vicinities, before being abruptly ejected 
away, so that exp erimental observations yield recurrent sequences of familiar patterns 
(JArtuso et all (fl990) : KersweU (:2005{ )). 



A yet unanswered issue concerns the initial conditions that are most suited to yield 
such unstable states. Traditional emphasis on so-called optimal perturbations may be 
misplaced. In fact, there is but a weak connection between the flow structures that 
grow the most during the linear transient phase and the chaotic flows found at large 
times. Such a connection for the case of the pipe flow concerns the so-called ed g e stat e 



which sits on a separatrix between laminar and turbulent flows ( Eckhardt et al\ ( 20071 ): 



Pringle fc Kerswelll ()2007l) ). This state, made up by two asymmetric vortices in the cros s- 



section, resembles the optimal disturbance of transient growth theory ( Bergstroml ( 1993[ )). 



For the motion in a square duct there seems to be no connection at all: the low-i?e turbu- 
lent flow, when averaged in time and space, is characterized by eight secondary vortices 
symmetric about diagonals and bisection lines. It seems reasonable to argue that such 
secondary structures represent the skeleton of the unstable periodic orbits, but the dis- 
turbances that grow the most in the linear transi e nt ph ase are formed by two vortices, 
symmetric about a diagonal (IGalletti fc Bottarol ( 2004 )). In both configurations, pipe 



and square-duct, the optimal perturbation is a stationary pseudo-mode, elongated in the 
streamwise direction, and not a travelling wave. This is a generic occurrence in wall- 
bounded shear fiows, and it does not bode well for the establishment of a simple, direct 
relation between small amplitude disturbances (excited in an initial receptivity phase) 
and finite amplitude wavelike states. Although nonlinear effects can be pinpointed right 
away as the culprit for the missing link between early stage of transition and late stages, 
there is scope for a receptivity analysis focussed on transiently growing initial condi- 
tions, followed by nonlinear simulations. A motivation for the search of wavelike struc- 
tures during t he early stages of transition is also provided by recent careful experiments 



(jPeixinho fc M ullin (2006)) on the reverse transition in pipe fiow, where modulated wave 



trains are found to emerge from long-term transients. 

The present paper starts by comparing the efficiency of optimal and sub-optimal per- 
turbations in triggering transition to turbulence at a value of the Reynolds number Re 
close to the threshold between laminar and turbulent fiow; it further shows that the 
turbulent motion oscillates around edge states which display an intriguing resemblance 
to optimal disturbances, before relaminarization occurs. Finally, an interpretation of the 
results is provided after projecting them onto a suitably defined phase space. 



2. Model configuration 

The incompressible fiow in a duct of square cross-section is an appealing configuration 
for the presence of geometrical symmetries capable to strongly constrain the patterns 
of motion. Countless studies have been devoted to the formation of secondary vortices 
in the turbulent regime (see iGavrilakia ( 1992 ) for a direct numerical simulation ap- 



proach) and, more recently, an attempt has been made to link the appearance of such 
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large-scale coherent states to the vortices appear i ng du r ing the initial, optim al transient 
phase of disturbance growth (JGalletti fc Bottarol (|2004[ ) : iBottaro et aLl(|2006h ). The lon- 



gitudinal la minar flow velocity compone nt has an analytic form U{y, z) available, for 
example, in lTatsumi fc Yoshimural ( 1990l ). with y and z cross-stream axes. After normal- 



izing distances with the channel height h, velocities with the friction velocity u^-, with 
u^ = ^-T-'jf-i time with /i/u,- and pressure with pu^, the following equations are found 
to govern the behaviour of the developed flow in an infinite duct: 

Ux+Vy+Wz = 0, 

ut + uux + vuy + wu^ = -px + -Rj^Aw -I- 4, 

Vt + UVx + VVy + WVz = -Py + T^^v, 
Wt + UWx + VWy + ww^ = -pz + -j^i-Aw, 

with A — dxx + dyy -t- dzz and Rct — u^ h/v. By using Ur as velocity scale we fix the 
pressure gradient, rather than the flow rate. 

An incompressible pseudo-spectral solver, based on Chebyshev collocation in y and 
z and Fourier transform along x, has been employed to solve these equations. For 
time-integration a third-order semi-implicit backward differentiation/ Adams-Bashforth 
scheme is used. In order to compute a pressure unpolluted by spurious modes, the pres- 
sure is approximated by polynomials (Pjv-2) of two units lower-order than for the velocity 
(P/v). Only one collocation grid is used, and no pressure boundary con dition is needed . 
The accuracy and stability properties of the method are discussed bv iBotellal (|l997l ). 



An adequate grid at Rer = 150 has been found to be composed by 51 x 51 Chebyshev 
points, with N^ = 128 streamwise grid points or 84 Fourier modes after de-aliasing. The 
a:-length of the domain has been chosen equal to An to accommodate a sufficiently large 
range of wavenumbers a, with periodic boundary conditions. A finer grid resolution has 
also been used for fully developed turbulent fiow with 71 x 71 x 256 physical grid points, 
or 71 X 71 X 170 spectral modes and a streamwise length L^ — GTip. The finer resolution 
run provides a slightly larger value of the threshold energy for transition, but integral 
quantities such as disturbance energy or skin friction factor are only marginally affected. 
Since the threshold value is, in any case, a function also of the shape of the initial con- 
dition, we do not deem indispensable to pursue expensive calculations to determine it 
exactly. For all cases studied, an adequate time step has been found to be Ai = 5 x 10"'*. 
The mean value of the generic function g is defined as g(y, z) — ^ ^^ / ^ g{x, y, z, t) dxdt. 
The addition of time averaging is necessary because of the finite (relatively low) stream- 
wise length. The bulk velocity is Ub = J u dydz and the centreline velocity is Uc = 
u(0.5,0.5). The friction factor for the square duct can be written as / = 8 u'^/U^. Some 
representative results are given in Table [U while the secondary flow field at Rcr = 150, 
averaged over forty units of time, is shown in figure [1] It displays a very regular pattern 
with eight vortices, despite the fact that no averaging over quadrants has been performed. 



3. Optimal perturbations 

The Navier-Stokes equations, linearized around the ideal laminar flow, are: 

iau + Vy + Wz = 0, 
Ut + iaUu -f vUy + wUz = —iap + {—a^u + Uyy + Uzz)/ReT, 
Vt + iallv = -py + {-a'^v + Vyy + Vzz)/ReT, 

f With this resolution we obtain an excellent match with the results bv lGavrilakij (|1992| ) for 
Re^ = 300. 
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f 


Uc/Ub Ret Rec 


laminar 0.018 


2.0963 3163 6630.4 


urbulent 0.0415 


1.53 2084 3188 



Table 1 . Comparison of some numerical values for laminar and fully-developed turbulent flow 
at Rer = 150. The subscrip t b ref e rs to bulk and c to centreline. The skin friction / given by 
the empirical correlation bv lJonea (| 19761 ) is / = 0.0481. 




Figure 1. Turbulent mean cross flow vortices and streamwise flow contours; isolines are 

spaced by 4 Ut. 



Wt + iaUw — —pz + (— a-w + Wyy + Wzz)/ReT, 

associated to boundary conditions u = v = w = Q on the walls; a is the streamwise 
wavenumber. The equations are integrated from a given initial condition at i = up to 
a final target time t = T. To identify the flow state at t = producing the largest distur- 
bance growth at any given T , a variational technique, based on the repeated numerical 
integration of direct a nd adjoint stability ecjuations, is used, coupled with transfer and 
optimality conditions ( Corbett fc Bottarol (2000)). The functional for which optimization 
is sought is based on an energy-like norm and reads: 



G{T) 



E{T) 



with E 



{u*u + v*v + w*w) dy dz, 



with * denoting complex conjugate. 

As a preliminary result the energy stability limit for this flow is determined; the min- 
imal Reynolds number below which disturbances decrease monotonically is found to be 
Rct = 23.22 for a = 0. In terms of Reynolds number based on centreline velocity and 
half channel height, this limiting value is equal to Rcc — 79.44. For comparison, in the 
case of Poiseuille flow, the Reynolds number is Rc r — 49, with a = and spanwise 
wavenumber /3 ~ 2 (c/. ISchmid fc HenningsonI (|200ll )). 

Then, we compute optimal perturbations at RCt = 150, for 3 different Fourier modes: 
a = 0, 1, 2. The gain and the corresponding cross-flow optimal disturbances at t = 
are presented in figure [2 A global optimal solution (G = 873.11) is found for a = 
at a time T — 1.31 and a comparable gain (G = 869.03) is found at a later time 
(T = 2.09) for a solution which is topologically different. These two solutions, which we 
call 'global optimals' represent streamwise vortices of vanishing streamwise disturbance 
velocity; they evolve downstream producing streaks of high and low longitudinal velocity. 
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Figure 2. Cross-flow velocity vectors for the optimal disturbances for a = 0,1 and 2. Note 
that for a 7^ the streamwise velocity fluctuations do not vanish at t = 0. 



A global optimal s t ate w ith two cells arranged along a duct diagonal, was obtained by 
Galletti fc Bottard (J2004l) in the context of a spatial, rather than temporal, optimization 
strategy. The existence of a four-cell optimal was not reported. It will be shown below that 
the two- and the four-cell states are very robust: when they are used as initial conditions 
for nonlinear simulations the flow trajectory is uncapable of evolving away from them 
towards a different topology which eventually exploits other symmetries, except when 
exceedingly large disturbance energies are given as input. 

For non-zero streamwise wavenumber, sub-optimal perturbations are found which take 
the form of modulated travelling waves. For a = 1, G = 381.41 at T = 0.384 and for 
a — 2, G — 237.05 at T = 0.246. To represent a wave train the temporal dependence of 
the generic disturbance can be written as the product of an exponential wave-part and an 
envelope function slowly modulated in time: q{x, y, z, t) = q{y, z, t)exp*"^^~'^*^ The phase 
velocity is found to be quasi-constant with time and close to the bulk velocity: c{a = 1) = 
1.1117 and c{a = 2) = 1.1536, both scaled with t/f,. The study of the temporal (rather 
than the spatio-temporal) evolution of disturbances is acceptable when flow structures 
travel at a well-defined speed within the duct. Such an approximation app ears to be 
reasonably well satisfied in experiments on equilibrium puffs in pipe flow (Hof e t ali 
(20051)), which are found to be advected downstream at speeds slightly larger than Ub, for 
values of the Reynolds number Ret — UhD/v [D pipe diameter) exceeding Rcb w 1800. 
Below such a threshold, turbulence can no longer be maintained autonomously. 



4. Nonlinear evolution 

Temporally evolving simulations have been conducted in a periodic duct of length 47r 
for a variety of initial conditions at Rcr = 150; such a val ue is close to the thre shold 
of self-sustained turbulence, as confirmed very recently by lUhlmann et al\ (2007). For 
each direct numerical simulation the initial state consists of the quasi-parabolic base flow 
profile plus a optimal or sub-optimal perturbation, normalized with prescribed energy Eq, 
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Figure 3. Skin friction for the four different optimal perturbations. For q = 0, the initial 
energy Eq is equal to 10"\ for a = l, Eo = 7.8 x 10"^ and for a = 2, Eo = 4.4 x 10"^. 



plus random noise. The complex amplitude of the noise in Fourier space varies between 
±10""'^''; the noise is necessary to fill the spectrum. 

In figure [3] the time evolution of the skin friction factor / is plotted for four initial 
conditions of different initial energies. When the simulations are initiated with one of the 
two global optimal solutions of figure [5] the ensuing behavior is uneventful (c/. figure [3]), 
and even for a rather large initial disturbance amplitude, £^o = 0.1, no instability appears 
to modiiy the fiow which returns slowly to the laminar condition with / — 0.018. The 
results are more suggestive when linear travelling waves are used to initiate the nonlinear 
computations. For the case in which the condition at i = is the sub-optimal state with 
a = 1 or 2, the initial growth of / is slower than for the global optimal case but, by time 
t = 2, a strong mean fiow deviation is created by nonlinear interactions, leading to a 
value of the friction factor which oscillates around the turbulent mean value / = 0.0415. 
It is notable that the energies Eq of the sub-optimal initial conditions sufficient to trigger 
transition are very much lower than 10~^. 

We have explored in more details the long time behaviour of the fiow when the sub- 
optimal initial perturbation with a = 1 is used to trigger transition; the results are 
condensed in figure |4l As a function of the value of Eq the Navier-Stokes calculations 
permit to recover the linear behaviour (when Eq ^ b x 1Q~^) or to depart from it. The 
threshold for transition is found for Eq — 7.8 x 10^^; below such an initial energy value 
the flow returns rapidly (within a few units of time) to the laminar state, above it the 
flow becomes turbulent. The energy gain of figure |3] shows that when Eq = 7.8 x 10~^ 
there is, at t = 0.8, a sudden increase of the fluctuation energy, likely linked to an insta- 
bility of the distorted mean fio w. There i s an in t eresting connection here with the recent 
theory of 'minimal defects ' Bottaro et all ( 20031 ): iBiau fc Bottard ( 20041 ): iGavarini et al. 
(|2004[ ): lBen-Dov fc CohenI (|2007| ). The secondary flow at t = 0.8 shown in the same figure 



displays symmetries about bisectors and diagonals, but this is not a generic occurrence 
and different initial conditions generate mean flow defects with other symmetries. For 
£^0 = 7.7 X lO^'^ the growth is followed by decay and rapid relaminarization (c/. also 
figure [5] left); when Eq — 7.8 x 10~^ the growth which starts at i = 0.8 is followed by 
a rapid filling of the spectrum with a peak in the intensity of fiuctuations at t — 1.8. 
Such a filling is made clear by figured] (right): first the modes with a = ti, n e N, 
grow because of quadratic interactions, then the modes with a = (2n — l)/2, n € N, 
emerge out of the random noise (visible in the figure after t = 5). Figure [5] (left) focusses 
on the modes a = and a = 2 which are first to be produced by nonlinearities. The 
streamwise-independent mode remains amplified for a long time because of lift-up effect, 
its appearance out of the a = ±1 fundamental mode is the analogue of the so-called 
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Figure 4. Evolution in time of the fluctuation energy, initial condition is given by the sub-op- 
timal perturbation with a = 1. Different curves correspond to different initial energy values: 
continuous line corresponds to i?o = 7.8 x 10"'^, dash-dotted line to E^ = 7.7 x fO~'^; dashed line 
to Bo = 5 X 10~^ and dotted line represents the linear computation. The energy is normalized 
with its initial value -Eo- The streamwise averaged vortices are drawn at times t — 0.8, 20, 50, 70 
and 80. 



oblique transition process in channel flow ( Schmid fc HenningsonI ( 19921 )). By time i = 3 
the turbulent flow can be considered as fully developed. Such an edge state persists un- 
til t « 78 (c/. [4]); it remains dynamically connected to the laminar base flow solution 
since relaminarization is abruptly reached at t « 80, after coalescence of the smaller 
vortices into a pair and then into a single large vortex. While turbulence is maintained, 
the secondary patterns displayed at i = 20, 50 and 70 in figure H] resemble the four-cell 
global optimal disturbance. At t — 20, two pairs of vortices are clearly visible in the 
cross-section; they are close to the two vertical walls whic h can thus be defined as 'ac- 
tive' since it is there that the turbulent wall cycle operates ( Uhlmann et all ( 20071 )). This 
4-cells state oscillates while maintaining remarkable coherence for some twenty units of 
time. Averaging over t yields the same flow pattern with four regular vortices discovered 
very recently bv lUhlmann et al\ (2007). Past t = 50, the active walls shift and the large 
scale vortices lean on the horizontal surfaces, despite the presence of smaller intermit- 
tent features that can be found near the left vertical wall at i = 50 and the right wall 
at t = 70. At i ~ 80 the flow relaminarizes; tests performed with different lengths of 
the computational domain show that the threshold Eq for transition remains unaffected, 
but relaminarization occurs later both for a shorter computational box (t « 90 when 
Lx = 27r) and a longer box (t « 150 when L^ = Btt), for the same numerical grid density. 
Interestingly, for the case of the short box the oscillations of the variables around the 
mean display enhanced amplitudes, highlighting the fact that a distorted dynamics could 
be caused by constraining too much the flow structures. 

For a geometrical description of the transition process, we choose the phase subspace 
spanned by two observables: the Reynolds number (based on bulk velocity) and the 
energy of the streamwise-averaged flow, noted Eu = 0.5 J ^U'^ + V^ + W'^dydz, with 
{U, V, W) the velocity vector after streamwise averaging. 

Figure |6] shows some trajectories between the laminar {Rcb = 3163, Ejj = 306.45) and 
the turbulent {Rcb = 2084, Ejj = 116) fixed points (the latter is simply defined by the 
temporal average of Ejj and Rcb in the fully developed turbulent regime) . The paths for 
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Figure 5. Left: temporal behavior of three different Fourier modes (a = 0, 1, 2) for the nonlin- 
ear simulation initiated by the optimal a = 1 initial condition with Eq = 7.8 x 10"'' (continuous 
lines) and Eq — 7.7 x 10""' (dashed lines), plus low amplitude noise. On the right the long 
time evolution of the fluctuations (a — 0.5, 1, 1.5, ...) for the case Eq = 7.8 x 10"'^ is repre- 
sented. The power density spectrum k„ is defined by k„ = -^ J ^(u*u)„ -I- (u*u)]v-n dydz, with 
u„ [oin, y, z, t) the a;-Fourier transform of u, from which the laminar profile has been substracted. 
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Figure 6. 'Phase diagram' representation of the transition process: time evolution in Reynolds 
number, based on bulk velocity, and mean flow energy Eij- The simulations start from the 
laminar flow solution (open circle in the left figure) plus a (sub-)optimal perturbation. The cross 
with arrows provides the approximate position of an unstable saddle node. The time averaged 
values of Eu and Kch in the fully developed turbulent regime are indicated by dashed lines in 
the zoom at right which refers only to the initial disturbance perturbations with a = 1. 



the initial conditions with a — Q correspond to homoclinic orbits in phase space. For the 
two cases which follow a non-trivial branch (a — \ and 2) the flow approaches the unstable 
saddle node (qualitatively sketched in the figure) before departing away from it along its 
unstable manifold. The trajectory in figure |6] (right) then starts circling around the point 
which characterizes the fully developed turbulent state, with orbits of increasing size. 
Before the end of the fifth orbit, it escapes through the unstable manifold of the saddle 
node towards the laminar fixed point. Each orbit has an ellipsoidal shape and is made up 
by two portions with very small local radii of curvature, around which the flow spends 
most of its time, and two long portions of large radii of curvature which are spanned very 
rapidly by the flow. Relaminarization is consistent with the emerging picture of shear flow 
turbulence which considers it as a tran sient event with a characteristic lifetime increasing 
exponentially with Reynolds number (|Hof et al\ (|2006l )). 
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5. Discussion and conclusions 

Although non- normality and transient growth are important issues, the traditional 
emphasis on optimal disturbances may be misleading when trying to predict the onset of 
transition. The results found here suggest that optimal initial perturbations in the form 
of steady streamwise vortices play mostly a passive role, while rapidly growing travelling 
wave packets have the potential to induce transition past a reasonably low threshold 
value of the initial disturbance energy. This suggests that an energy-like functional is 
possibly not the most pertinent objective of the optimization procedure. 

In the simulations performed here we have found that key to transition is the estab- 
lishment of a disturbed mean flow profile (mode a — 0), susceptible to an exponential 
or algebraic instability which causes enhanced growth of the travelling wave mode with 
a = ±1. Such an early stage can be described by a weakly nonlinear triadic interaction 
model. 

Once turbulence is established, its sustainment depends on non-normality, capable 
to produce rapid transient growth of the disturbance energy, and non-linearity, which 
enables directional redistribution of the amplified disturbances. The coupling between 
transient growth and directional redistribution of energy causes the flow to orbit around 
the turbulent fixed point in the phase-space of figure [HI Near the edge of chaos the 
lifetime of turbulence is finite and beyond a given value of t (function of the Reynolds 
number and of the streamwise dimension of the periodic domain) relaminarization occurs. 
Interestingly, the edge state resembles t he op t imal disturbance; this fact has recently been 
observed in pipe flow ( Eckhardt et al\ ( 2007 ): |Pringle fc Kerswelll ( 2007 )) and still awaits 



for an explanation. 
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